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Abstract. We study the evolution of embedded protoplanets in a protostellar disk using very high resolution 
nested-grid computations. This method allows us to perform global simulations of planets orbiting in disks and, 
at the same time, to resolve in detail the dynamics of the flow inside the Roche lobe of the planet. The primary 
interest of this work lies in the analysis of the gravitational torque balance acting on the planet. For this purpose 
we study planets of different masses, ranging from one Earth-mass up to one Jupiter-mass, assuming typical 
parameters of the protostellar disk. The high resolution of the method allows a precise determination of the mass 
flow onto the planet and the resulting torques. The obtained migration time scales are in the range from few 
times 10^ years, for intermediate mass planets, to 10® years, for very low and high mass planets. Typical growth 
time scales depend strongly on the planetary mass, ranging from a few hundred years, in the case of Earth-type 
planets, to several ten thousand years, in the case of Jupiter-type planets. 
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1. Introduction 

During the past five years radial velocity studies 
have allowed the detection of planetary companions 
around other main-sequence stars. Until now about 
sixty so-called “extrasolar planets” have been discov¬ 
ered, which orbit their stars within a distance of a 
few AU. A recent catalog of extrasolar planets, in¬ 
cluding their orbital characteristics, is provided by 
Butler et al. ( pOOl ) and up-to-date versions can be 
found at http://www.obspm.fr/encycl/encycl.html 
and http://exoplanets.org/, maintained by Jean 
Schneider and the Department of Astronomy at UC 
Berkeley, respectively. 

In contrast to the solar system, these new planets dis¬ 
play quite different orbital properties that challenge the 
accepted formation scenario for solar planets. The ma¬ 
jor differences are their high minimum masses (up to 17 
Jupiter-masses), their proximity to the central star (a frac¬ 
tion of the Sun-Mercury distance) and their high eccen¬ 
tricities (up to 0.7). 

One of the main problems to deal with is the very 
close distance of massive planets to their parent star. The 
formation of Jupiter-type planets at these locations is, on 
theoretical grounds, very unlikely. First of all, from purely 
geometrical arguments, the matter reservoir of the sur¬ 
rounding disk is too little so that a planet could hardly 


accrete its mass. Second, the temperatures within the disk 
are too high for a rocky core to condense easily. 

For these reasons it is generally believed that planets 
have formed from disk material further out, at distances of 
several AU from the star, and have then migrated to their 
present positions. This radial motion of the planet through 
the disk is primarily caused by gravitational torques acting 
on the planet. The presence of the planet in the disk dis¬ 
turbs the disk gravitationally, creating spiral density wave 
perturbations, which emanate from the planet through the 
disk. Hence, the disk is no longer axisymmetric which re¬ 
sults in a net torque on the planet. The sign and magni¬ 
tude of the vertical component of the torque determines 
the direction and efficiency of the radial migration. 

While initial fully non-linear hydrodynamical numer¬ 
ical computations of embedded pla nets a ssumed a fixed 
circular orbit of the planet (Kley 1999; Bryden et al. 
1999; Lubow et al. 1999), more recent simulations took 
into account the back reaction of the disk and allowed for 
a change in the parameters of the planetary orbit (Kley 
^000 ; Nelson et al. 2000). For a Jupiter-mass planet and 
typical parameter values for the disk, the obtained orbital 
decay time is about 10® years, which agrees reasonably 
well with previous estimates based on analytic linear the¬ 
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ories (Goldreich & Tremaine 1980; Ward 1997). 

The majority of the computations, performed so far, 
have used a single grid which resolves the Roche lobe 
of a Jupiter-mass planet only with very few grid cells. 
Recently, Cieciel§g et al. (2000a, 2000b|) used an Adaptive 
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Mesh Refinement method to resolve the immediate sur¬ 
roundings of the planet, but they didn’t give any estimate 
of the mass accretion rate and magnitude of the gravi¬ 
tational torque. On the other hand, Armitage (2001) re¬ 
duced the overall simulated region achieving a better res¬ 
olution. However, also in this case the Roche lobe is only 
resolved by a few grid cells because of the low mass of the 
investigated planet. 

In this paper we aim at the structure and dynamics 
of the gas flow in the close vicinity of the planet, while 
performing global disk simulations. In order to obtain the 
necessary high spatial and temporal resolution, we use a 
nested-grid formalism which allows an accurate compu¬ 
tation of the mass flow onto the planet and the acting 
torques. 

In the next section we layout the physical model fol¬ 
lowed by a description of the numerical method (Sect. 

We describe the setup of the various numerical models in 
Sect. The main results are presented in Sect. ^ and our 
conclusions are given in Sect. 


2. Physical model 


For the purpose of this study, we assume that the opening 
angle of the protostellar accretion disk is very small. We 
describe the disk structure by means of a two-dimensional, 
infinitesimal thin, model using vertically averaged quanti¬ 
ties, such as the surface mass density 



where p is the regular density. We work in a cylindrical co¬ 
ordinate system (r, z) whose origin is fixed at the center 
of mass of the star and the planet, and where the plane of 
the disk coincides with the z = 0 plane. 

The gas in the disk is non-self-gravitating and is orbit¬ 
ing a protostar having a mass M*. = 1 Mq . The total mass 
of the disk Md, within the simulated region, which extends 
from 2.08 to 13 AU, is 3.5 x 10“^ Mq. Embedded in this 
disk there is a massive protoplanet with a mass Mp, which 
ranges from one Earth-mass (Mg) to one Jupiter-mass 
(Mtj,), depending on the considered model. The planet is 
assumed to be on a fixed circular orbit throughout the 
evolution. We employ a rotating coordinate system, coro¬ 
tating with the planet, whose azimuthal position is kept 
constant at (/Jp = tt. The angular velocity Q of the rotating 
frame is then given by 


U — Up — 


G (Mi, + Mp) 


( 1 ) 


where G is the gravitational constant and a is the semi¬ 
major axis of the planet’s orbit. 

The evolution of the disk is given by the two- 
dimensional (r, ip) continuity equation for S and the 
Navier-Stokes equations for each of the two components of 
the velocity field u = {ur,u^). Thus, the set of equations 
read 

—+V.(S«) = 0, (2) 


+ V • (S m) = E r (w +11)2 -- E ^ (3) 


9[Er2 (a; -|- H)] 

di 


o dP 

+V-[Er2 {cu+n)u]=- — -E—+U.{4) 


Here ui = u^pj r is the angular velocity and P is the ver¬ 
tically integrated (two-dimensional) pressure. The gravi¬ 
tational potential d), generated by the protostar and the 
planet, is given by 


^' = $* + d>p 


GM* GMp 


( 5 ) 


where and Tp are the radius vectors to the star and the 
planet, respectively. The effects of viscosity are contained 
in the terms fr and which give the viscous force per unit 
area acting in the radial and azimuthal {fp/r) direction: 


1 d{r Srr) 1 dSrp Spp 
r dr r dip r ’ 


_ 1 9(r2 Srp) dSpp 
^ r dr dip 


Since we assume a zero bulk viscosity C, a constant kine¬ 
matic viscosity v and we do not include any artificial vis¬ 
cosity in our numerical models, the relevant non-zero com¬ 
ponents of the three-dimensional viscous stress tensor S 
are 


dur 1 
~ 3 


Spp = 2vT, 


Srp — E 


dip r 

1 dur 

- -J, -k f 

r dip 


V-ttj , 

(6) 


( 7 ) 

duj\ 

(8) 

dr I 


where the divergence of the velocity field can be written 
as 

V 

r or dip 

A more general form of the relations §), ( 0 ) and (||), 
within the two-dimensional cylindrical approximation, is 
given for example in Kley (1999). 

In the set of equations above we have omitted the en¬ 
ergy equation because in this study we will be concerned 
only with a relatively simple equation of state which does 
not require the solution of an energy equation. We shall 
use an isothermal equation of state where the surface pres¬ 
sure P is related to the density E through 

P = c2e. (9) 


The local isothermal sound speed Cg is given here by 

Cs — ^Kepi 

r 

where x;Kep = \/G M^fr denotes the Keplerian orbital 
velocity of the unperturbed disk. Equation (0) follows 
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from vertical hydrostatic equilibrium. The ratio h of the 
vertical height H to the radial distance r is taken as a 
fixed input parameter. With this choice, the temperature 
T is proportional to 1/r. Here we use a standard value 

h = — = 0.05, 
r 

which is typical for protostellar accretion disks having a 
mass inflow rate of M « 10“^ Mq yr“^. With this value of 
h, our kinematic viscosity coefficient is equivalent to a = 
4 X 10“^ at the radial position of the planet {v = a Cg H). 

Since the mass of the planet is very small in comparison 
to the mass of the star, because we always use here a 
ratio q = Mp/ M* smaller than 10“^, the center of mass 
is located very close to the position of the star. In the 
following we will often identify the radial distance from 
the origin of the coordinate system with the distance from 
the central star. 



Fig. 1. Face-up projection of a three-level grid system in 
cylindrical coordinates. On the finest subgrid {I = 3) the 
linear spatial resolution is four times as large as it is on 
the main grid {I = 1). 


3. Numerical method 

In order to study the planet-disk interaction, we utilize a 
finite difference method to solve the hydrodynamic equa¬ 
tions outlined in the previous section. As we intend to 
achieve a very high resolution around the planet, we use a 
nested-grid technique. Both requirements are provided by 
an earl y FO RTRAN-Version of Nirvana (Ziegler 1997 ; 
Ziegler 1998 ), which is a 3D, nested-grid, MHD code, 
based on a covariant Eulerian formalism. The relevant 
equations are solved on a mesh structure having a con¬ 
stant spacing in each direction. For the present purposes, 
the code is used in a pure hydrodynamic mode, adopt¬ 
ing a cylindrical reference frame where the z-dimension is 
switched off. 

Nirvana uses a spatially, second-order accurate, ex¬ 
plicit method. The advection is computed by means of the 
second order monotonic transport algorithm, introduced 
by van Leer (1977), which guarantees global conservation 
of mass and angular momentum. It is first-order accurate 
in time. The viscosity part was added and is treated ex- 
plicitely. The stress tensor has been implemented for a 
Newtonian fluid according to the Stokes hypothesis, i.e. 
with a bulk viscosity C = 0. 

This code has been already empl oyed i n similar com- 
putat ions, either in 2D (Nelson et al. 2000 ) and 3D (Kley 
et al. 2001 ), but always in a single-grid mode. 


A similar numerical scheme has been adopted, for as- 
trophysical simulations, by a number of authors. Ruffert 


(1992) used this approach to investigate the collision be¬ 
tween a white dwarf and a main sequence star. In his paper 
the numerical method is explained in detail. Yorke et al. 


(1993) and Burkert & Bodenheimer (1993) simulated the 


collapse of a protostellar cloud. An application to flux- 
limited radia tion hydrodynamics can be found in Yorke & 
Kaisig ( 1995|) . 

The method relies on the basic idea that, whenever a 
greater resolution is needed in a designated region, a finer 
subgrid is located inside the main grid (the one covering 
the whole computational domain). If the resolution is not 
high enough yet, another subgrid may be placed on the 
underlying one. Since any subgrid can host a finer sub¬ 
grid structure, a grid hierarchy is generated, also called 
“system of nested grids”. In principle there is no limit to 
the degree of nesting. A three-level hierarchy is shown in 
Fig. 0. 

The necessary equations are then integrated, indepen¬ 
dently, on every grid level. However, two neighbor sub¬ 
grids must exchange the necessary information whenever 
the integration proceeds from one grid level to another. 
Restrictions are imposed on the time step only because, 
for stability reasons, the Courant-Friedrichs-Lewy (CFL) 
condition must be fulfilled during each integration, on each 
level. 


3.1. Nested-grid technique 

This technique is particularly useful when very high local 
resolution is required at specific and predefined points of 
the computational domain. In our situation, this allows 
us to simulate both the overall behavior of the disk and 
the immediate surroundings of the planet. Since this kind 
of numerical approach is quite new for the calculation of 
disk-planet interaction, we describe the method to some 
extent, but only referring to our particular and specific 
case. 


3.1.1. Basic integration cycle 

In our calculations we use the smallest possible refinement 
ratio: 

Ar{l + 1) Aipjl + 1) ^ 

Ar{l) A(p{l) ’ ^ ’ 

where Ar(l) and Aip(l) represent the mesh discretization, 
along each direction, on the grid level I (here I = 1 identi¬ 
fies the main grid). In order to analyze a complete integra¬ 
tion cycle, let’s suppose we have a three-level hierarchy at 
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an evolutionary time t. The cycle will be completed when 
on each grid the system has evolved for the same time: 


1. We always start the integration from the finest level. 
During the first step of the cycle, this grid is evolved 
for a time interval 


Ati(3) = min 


At^ 




where Atf^^{l) represents the time step resulting from 
the CFL criterion applied to the level I, after its latest 
integration. Thus, Ati(3) accounts for the CFL stabil¬ 
ity criterion on the whole set of grids. 

2. The third grid has to be integrated once more because 
of Eq. &■ Since after the first step this was the only 
level to evolve, we only have to check the new CFL 
time step for this grid, Af2^'^(3). Then it can move 
further in time for an interval 


At2(3) = min 




3. Now level 2 can be integrated for a time 
At3(2) = Ati(3) + At2(3) < Aff’"L(2), 


so that numerical stability is automatically assured. At 
this point of the cycle the first information exchange 
takes place: the solution on the level 2 is corrected via 
the more accurate solution of the level 3; the boundary 
values of the level 3 are updated by using the solution 
of the level 2, which covers a larger domain. These 
fundamental interactions will be described later. 


We have just seen that a level I -I-1 has to be visited two 
times as often as the level 1. Then the next two cycle steps 
will be similar to the first two, provided that the appro¬ 
priate CFL time steps are employed to compute Af4(3) 
and Af5(3). During the sixth step, the level 2 evolves for 
At6(2) = At4(3)-t-Af5(3). Eventually, during the seventh 
step, the integration of the level 1 is performed using a 
time step Af 7 (l) = At 3 ( 2 )-f At6(2) which is, by construc¬ 
tion, smaller than Aff^'"(l). The cycle is now complete 
and each grid level has evolved for the same amount of 
time At 7 (l). The entire cycle sequence is schematically 
sketched in Fig. ||. 

In general, within this kind of cycle, a level I is inte¬ 
grated 2^“^ times. 


3.1.2. Downward information transfer 

We already mentioned that after the third step of the it¬ 
eration cycle, grids 3 and 2 have to exchange some infor¬ 
mation. In general, this exchange must occur every time 
the grid level I evolves to the same time as the level Z -|- 1. 
Because of the higher resolution, we assume the solution 
of the level Z -I-1 to be more accurate than that of the level 
Z. Therefore, the fine-grid solution replaces the coarse on 
the common computational domain. Whatever the level in 
the hierarchy is, the frame formed by the first and the last 



Fig. 2. Scheme of a complete integration cycle for a three- 
level grid system. Arrows indicate the direction of infor¬ 
mation transfer when integration proceeds from a level to 
the next lower one. Straight arrows stand for the solution 
updating process on the levels 2 and 1. Bow-arrows indi¬ 
cate the data transfer for setting the boundary quantities 
on the levels 3 and 2. 



Fig. 3. Interface between a subgrid and its host. The light- 
colored zone marks those cells containing boundary values 
needed for the subgrid integration, the “ghost cells”. The 
darker region refers to the so-called “active zone”, where 
values are effectively computed on the subgrid. The thick 
line, which separates the previous regions, encloses the 
grid cells on the coarse grid whose content is replaced by 
the more accurate one coming from the subgrid. 


two grid cells are ghost cells (see Fig. This indicates 
that they contain the boundary values necessary to per¬ 
form the algorithm integration. Ghost cells of level Z -|- 1 
do not contribute to the updating process of the solution 
of level Z. 

The replacement procedure is straightforward: a sur¬ 
face weighted average, using the nearest fine values, substi- 
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Fig. 4. Surface weighted average of the surface density. 
A coarse cell is shown along with the four fine cells it 
comprises. As a scalar, the surface density is cell-centered 
within the appropriate refinement (dots). represents 
the new, interpolated, value of the coarse cell; Ef are the 
fine interpolating quantities. Because of the fixed value of 
Aip{l), Al = A 4 and A 2 = A 3 . 


tutes the coarse quantity. For example, referring to Fig. 
the averaged coarse density (E®) is 


.c ^ E.SfA, ^ (EF + Ef)Ai + (Ef + Ef)A2 
y~! .- Ai 2 (Al -i- A2) 


In our case, since the advected quantities are the linear ra¬ 
dial momentum density E Ur , and the angular momentum 
density E cj, these are the interpolated quantities, along 
with E. Since velocities are centered at the sides of a cell 
(see Figs. ^ and this average is a little more complex 
than the previous one and requires six terms. Indicating 
with the coarse value of the linear momentum density 
to be interpolated and with uj the surrounding fine grid 
values, we have: 


_ (Mi + ^6 ) ^1 + (’^2 + ^5 ) ^2 + (^3 + ^4 ) A 3 


2 (Al -|- A 2 + A 3 ) 


In a successive step, velocities are retrieved. No other 
quantity needs to be interpolated because the potential is 
fixed and no energy equation is solved. 

In order to guarantee global mass and momen¬ 
tum conservation for the whole hierarchy, we have 
to make sure that the momentum flux components 
(Ar Huj./At, r Aip E w/ At ), across the border between 
a subgrid and its host grid (indicated in Fig. ||), are equal 
in both level solutions whenever the grid has evolved for 
the same time as the subgrid. However, each grid evolves 
independently and for a time interval different from that 
of the lower one. Thus, even after the solution updat¬ 
ing process described above, the amount of momentum 
flowed across the borders might not coincide in the re¬ 
spective solutions. To remove this possible discrepancy, 
at the coarse-fine grid border, these quantities are taken 
from the fine grid integration. In Fig. the situation for 
the azimuthal momentum flux is depicted. Two fine cells 
participate in this process. Referring to the integration 



Fig. 5. Surface weighted average of the radial momen¬ 
tum density u = Eu^. Two coarse cells are drawn (thick 
lines); Fine cells are delimited by thin lines. Because of the 
staggered structure of the grid, vectors are face-centered 
within the cell. 


cycle traced in Sect. 3.1.1, fl{l) represents the value of 
the quantity r Aip E ui/At, at the grid-grid interface lo¬ 
cation, as computed during the fc-th cycle step on level 1 . 
An additional index (j = 1 , 2) is needed to identify the ra¬ 
dial position of the two fine cells involved (for example, on 
level 3), but it does not concern the coarse grid quantity 
to be replaced (on level 2 ). 

Suppose we are at the end of the third cycle, when 
the first interaction, between levels 2 and 3, occurs (first 
straight arrow in Fig. |). Because of the refinement ratio 
established by Eq. (p), quantity / 3 ( 2 ) will be reset as: 


/3(2) = 2 


1 Ati(3) E, /f(3) + At2(3)E,/l(3) 


Ah{2) 


(14) 


This correction is accounted for directly while performing 
the advection of radial and angular momenta. 


3.1.3. Upward information transfer 

The boundary conditions on the main grid are usually 
imposed depending on the physics and geometry of the 
problem: symmetry, periodicity etc. In the case of a sub¬ 
grid, boundary values must be attached, in some way, to 
the values of the underlying grid. This point turned out 
to be extremely delicate for our calculations. In Nirvana, 
the boundary values of a certain level I were just set by 
means of a linear interpolation of the quantities on the 
lower level I — 1. Because of the strong variations in den¬ 
sity and velocity, due to the formation of shock fronts in 
our simulations, this method fails and produces numerical 
inconsistencies. 

Therefore, we raised the order of the interpolation. 
However, this introduces another potential trap. In fact, a 
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Fig. 6. Momentum flux correction scheme. Momentum 
flux components are centered as velocity field compo¬ 
nents, whose locations are shown in the upper left corner. 
/(2) Ar(2) = /a(2) At3(2) Ar(2) represents the amount of 
angular momentum, flowed across the coarse cell border, 
during the third cycle step (see Sect. ^.l.l| ). /^(3) Ar(3) = 
/j (3) Ati(3) Ar(3) -I- (3) At2(3) Ar(3) represents the 

same quantity transported, during the first two cycle 
steps, across the j = 1 fine cell border. The coarse quan¬ 


tity is replaced by 


f(3) + /2(3) 


Ar(3). 


high-order interpolation (higher than the first order) is not 
monotonic and can produce a new minimum. This is not 
acceptable since, for example, the density is a non-negative 
quantity. Then the interpolating function should be mono- 
tonised. In order to handle this problem we h ave used 
the same approach as described in Ruffert ( 1992| ), that is 
by employing the monotonised harmonic mean (van Leer 


1977). 


Then, if we have a function sampled at a; — Arc, x and 
X -\- Ax, with values gi,, g and g^, respectively (as shown 
in Fig. 1^), the averaged value at x -|- e is 

{g - 5 l ) ( 5 r - g) 


2 e 

g<^ = g + ^ max 


ffR - 5 l 


-,o 


(15) 



Fig. 7. Behavior of an harmonic mean against a geomet¬ 
ric mean. Grey filled circles indicate the harmonic values 
assumed by g^. Open circles indicate the corresponding 
values when the geometric mean is applied. The plot in¬ 
tends to simulate a shock front. If we keep the left value 
unchanged and make the right one g^ four times as large, 
the geometric mean value will increase proportionally. In 
turn, the harmonic mean will change so slightly that its 
variation is not visible on this plot. 


lies on a coarse cell border or whether it does not. In 
the first case, either Sr or e^p is zero. Then only three 
coarse values participate in the average, along either the 
azimuthal or the radial direction. In the second case the 
interpolation proceeds exactly as explained for the surface 
density. 


4. General model design 


The main goal of this study is the investigation of the 
characteristic features of the hydrodynamic flow within 
the Roche lobe of the planet. This means that we must be 
able to resolve a characteristic length, the Hill radius: 





(16) 


provided that —Ax/2 < e < Ax/2. 

If we adopt this kind of average on a 2D-mesh, each 
interpolation generally involves 3x3 coarse quantities. It 
proceeds by averaging the selected coarse values along a 
certain direction, three at a time. This results in three new 
quantities. A further harmonic average of these, along the 
other direction, generates the subgrid boundary value at 
the correct position. Fig. ^ shows how the procedure works 
in the case of the surface density. 

As density is cell-centered, er and are always one 
fourth of the coarse grid linear size. Since the averaging 
process is performed for each direction separately, it is not 
affected by the metric of the mesh, i.e. it is performed as 
if a Cartesian grid were used. 

For the velocity components we have to distinguish two 
different cases (see Fig. H): whether the boundary value 


Moreover, we intend to do that for a variety of planet- 
to-star mass ratios. In order to reach such resolutions, we 
build a series of grid systems having the planet located 
approximately in the middle, at each level of the hierarchy. 
Smaller planets require higher degrees of nesting. 

From now on, we refer to non-dimensional units. All 
the lengths are expressed in units of the semi-major axis 
a. This is constant because the planet moves on a fixed 
circular orbit, with radius 


Masses are in units of the central stellar mass and time 
is given in units of the planet’s orbital period. However, 
in order to convert them into conventional physical units, 
we assume that a = 5.2 AU and, as already mentioned, 
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Fig. 8. Harmonic average: surface density boundary val¬ 
ues. The light-colored zone indicates subgrid ghost cells. 
Darker region belongs to the active subgrid zone. Nine 
coarse values (big black circles) are engaged in the av¬ 
erage. The value to be interpolated is shown as a small 
black circle. During a first step, Eq. ( p^ is applied three 
times along the (/^-direction. Three new values are gen¬ 
erated (grey small circles), having the correct, final az¬ 
imuthal coordinate. These are averaged to obtain the final 
value at the correct radial location. Small open circles re¬ 
fer to the other boundary quantities whose value depends 
on the same coarse quantities. 


M* = 1 Mq. This implies that one planet orbit takes 11.8 
years. 

The whole azimuthal range of the disk is taken into ac¬ 
count by considering a computational domain represented 
by 2 7r X [rin,rout], where Hn = 0.4 and rout = 2.5. This is 
covered by a 142 x 422 mesh (main grid), allowing a res¬ 
olution such that Ar(l) = A(/5(l) = 0.015 and constrain¬ 
ing the resolution on each other grid level, according to 
Eq. (0). The size of any higher grid level is 64 x 64. 


4.1. Smoothing of the potential 


The perturbing action of the planet is exerted via its grav¬ 
itational potential 4>p. From a numerical point of view, it 
is usually smoothed in order to prevent numerical prob¬ 
lems near the planet. Thus, we re-write the denominator 
of <I)p in Eq. (i as However, the smooth¬ 

ing length 5 cannot be the same all over the grid system 
because of the different grid size involved at each level. 
Then we use the following grid-dependent length: 


S{1) = min 



(18) 


where 

X{1) = ^Arm)+rlAip^l) ~ V2Ar{l). 


(19) 


The value of constant part, in Eq. 


has become a kind 


of standard in single grid simulations (Kley 1999; Lubow 


Fig. 9. Harmonic average of radial velocity component. 
Two cases are shown. In one case the position of the ve¬ 
locity component lies on a coarse grid border (thick lines), 
requiring only three coarse quantities (one single average). 
In the other case nine coarse quantities are needed as for 
the density interpolation. The only difference is that Sr is 
one half of the radial coarse size. 


et al. 1999; Kley 2000). This is always used on the main 
grid whereas the grid dependent part always prevails on 
the highest grid level. This choice results in a very deep 
potential in the immediate vicinity of the planet. 


4.2. Accretion onto the planet 

The presence of the planet affects the nearby disk den¬ 
sity also because it can accrete matter. Planet accretion 
is accounted for by removing some mass from the region 
defined^ by |r — rp| < Kac- Since mass is removed from 
the system after each integration step, the evacuation rate 
depends on an input parameter /tev as well as on the in¬ 
tegration time step. The details of the accretion process 
are given in Fig. This removal is accomplished only on 
the highest (finest) hierarchy level and the removed mass 
is not added to the dynamical mass of the planet, but 
just monitored. A standard value Kev = 5 is used, while 
we set Kac between 8.0 x 10“^ and 9.4 x 10“^ i?H- These 
values are such that only few grid cells, on each side of 
the planet, are involved in the gas accretion, making it 
a locally confined process. For the Jupiter-mass planet, 
~ 12 X 12 cells participate. The lowest number of em¬ 
ployed cells is ~ 8x8, which is used for the smallest planet 
(Mp = 1 Mj). The largest is ~ 18 x 18 and is adopted for 
a planet with Mp = 0.5 (160 Afj). This circumstance 

is due to the high numerical resolution of the model (see 
Table I). 

^ Using this notation, we refer not only to the radial extent 
of the accretion region but also to the fact that the region is 
centered at rp. 
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Fig. 10. Planet mass accretion scheme. Wherever a cell 
center lies in the dark-colored zone, after each integration 
its density is lowered by an amount AE such that AS/S = 
2 TT Kev At, where At is the integration time step interval. 
If the cell center falls into the light-colored zone AS is 
only one third of that value. For a typical main grid time 
interval and Kev = 5, we would get AS/S « 1.2 x 10“^. 
Anyway, this evacuation process is performed only on the 
highest grid level for better accuracy. 


Note that the sphere of influence of the accretion pro¬ 
cess consists of a region typically 10“^ i?H which is smaller 
than the radius of the protoplanet, which is assumed to ftll 
its Roc he lob e during its growth p hase (Bodenheimer & 
Pollack 1986 ; Tajima & Nakagawa 1997| ). As the present 
study does not take the energy equation into account, it 
precludes a detailed treatment of the intern al structure of 
the protoplanet (see Wuchterl et al. 2000| and references 
therein, and also Fig. below), hence the inferred accre¬ 
tion rates may still be unreliable. 


4.3. Initial and boundary conditions 

The initial density distribution is proportional to 
However, we superimpose to this an axisymmetric gap 
around the planet, obtained by an approximate balance 
of the viscous torque and the gravitational torque due 
to the planet (P. Artymowicz, private communication). 
Figure ^ shows the surface density, at t = 0, for some 
selected planet masses. The initial velocity field is that of 
a Keplerian disk. 

As boundary conditions, periodicity is imposed at (/? = 
0 and ip = 2 TT. We allow matter to flow out of the compu¬ 
tational domain at the inner radial border (ri„) whereas 
we set reflective boundary conditions at the outer radial 
border (rout)- The angular velocity is set equal to the un¬ 
perturbed Keplerian value ^Kep = y/GM.*./r^, both at 
Tin and rout- 



Fig. 11. Initial surface density distribution for q = 10 
q = 10“"^ and q = 10“®. S(r, ip) = S(r) at t = 0. The gap 
greatly reduces when q gets small enough. 


For low-mass planets (Mp < 10 M^), boundary condi¬ 
tions should not affect much the system evolution because 
density waves damp before reaching r = rout and are very 
weak when they reach r = rin. For more massive planets, 
some reflection is seen at the radial outer border. Further, 
the torque exerted by the planet pushes the inside-orbit 
material inwards. As a result, because of the open bound¬ 
ary at r = rin, the inner disk is partially cleared. The 
higher the planet mass is the stronger these effects ap¬ 
pear. 


4.4. Model specifications 

In this paper, we are mainly interested in investigating 
how disk-planet interactions change by varying q, the 
planet-to-star mass ratio. Table ^ summarizes the values 
of q used in different models, along with the number of 
grid levels employed. For reference, a measure of the lin¬ 
ear resolution in units of the Hill radius is given 
for the highest level as well. Some models have different 
prescriptions than the ones outlined above, as specified in 
the Table |^. 

Few models may deserve some comments. Model 
Elen2 and WproI aim at checking whether results from 
CiROl and Ciro2 (respectively) are resolution-dependent. 
This test is negative, as we show in the next sections. Since 
the planet position (rp, ipp) can fall anywhere within a grid 
cell (according to the value of q in Eq. 0 and to the def¬ 
inition of the grid), some asymmetries could arise. These 
might have some effects on the finest levels, due to the 
small value of the smoothing factor S. In order to achieve 
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Table 1. Model specific parameters: q = Mp/M*; ng is 
the number of grid levels; ^ = Ar(ng)/Rn is the normal¬ 
ized grid resolution on the finest level. Each model is let 
evolve, at least, till 200 orbits. 


Model 


t 

f 

ng 


C 

Notes 

CiROl 

1.0 

X 

10"^ 

5 

1.3 

X 10"^ 


Ciro2 

1.0 

X 

10"^ 

6 

1.5 

X 10"^ 


Ciro3 

1.0 

X 

10"® 

7 

1.6 

X 10"^ 


PeppI 

3.0 

X 

10"® 

7 

2.3 

X 10"^ 


Pepp2 

1.5 

X 

10"® 

7 

1.4 

X 10"^ 


Pepp3 

3.0 

X 

10"® 

7 

1.1 

X 10"^ 


Pepp4 

6.0 

X 

10"® 

7 

8.6 

X 10"® 


WproI 

1.0 

X 

10"^ 

7 

7.3 

X 10"® 


Wpro2 

1.0 

X 

io-'‘ 

6 

1.5 

X 10"^ 

(a) 

GinoI 

1.0 

X 

10"® 

7 

1.6 

X 10"^ 

(6) 

Gino2 

2.0 

X 

10~'‘ 

6 

1.2 

X 10"^ 


Gino3 

5.0 

X 

10"'‘ 

6 

8.5 

X 10"® 


ElenI 

1.0 

X 

10"® 

5 

1.3 

X 10"^ 

(=) 

Elen2 

1.0 

X 

10"® 

6 

6.8 

1 

O 

X 



Same as Ciro 2, but Kev = 10; 


Same as CiROS, but the planet is symmetrically 
placed with respect to the grid cells. 

Same as CiROl, but reflective boundary conditions 
are set at r = Hn. 


a complete symmetry, in the model GiNOl the planet is 
placed at the corner of a main grid cell (i.e. the intersect¬ 
ing point of four grid cells). This property is such that, 
on every other gird level, the planet always sits on the 
cross-point of four grid cells. CiROS and its counterpart, 
GinoI, give almost identical results. Since in the Jupiter- 
mass case the inner-disk is greatly depleted, model ElenI 
was run to evaluate the influence of its presence on the 
gravitational torque. For this reason, in such model we 
prevent matter from draining out of the inner radial bor¬ 
der by setting reflective boundary conditions. 


5. Results 

Hereafter we mainly discuss three models, namely CiROl, 
Ciro2 and GiR03. We concentrate on them because 
they cover a mass range from 1 down to 3.2 Mg. 
Nevertheless, whenever required by our discussion, we 
mention other particular models. Some results, concerning 
the whole set of models given in Table |^, are presented as 
well. 


5.1. Overall flow structure 


Large-scale interactions (whose effects extend over 2 tt in 
azimuth and cover a large radial extent) of a Jupiter- 
mass planet with the surrounding environment, have al¬ 
ready been t reated num erical ly in a numbe r of p apers 
(Arty mowicz 1992 Kley 1999 ; Bryden et al. 1999| ; Kley 
200 1|). A n example of large-scale features can be seen in 
Fig. where the comprehensive result of a nested-grid 
computation is displayed. 


Planets with a lower value of q should perturb the disk 
less and have a weaker large-scale impact on it. For these 
reasons we discuss only the medium ( 1 ^ — Tpl ~ i^n) and 
small (|i— Tpl <C Rh) scale effects of such interactions. 

Nonetheless, it’s worthwhile noticing that large-scale 
structures are clearly visible in our smallest mass models. 
In Giro3 (Mp = 3.2 Mg), a trailing density wave, ema¬ 
nating from the planet, spirals around the star for about 
Jtt, vanishing approximately at r = 2. In PeppI (Mp = 1 
Mg), a similar feature spirals for almost 2 7 r, disappearing 
at r = 1.5. Although we did not investigate this issue any 
deeper, it may happen that results from local simulations 
could be influenced by not taking into account entirely 
such global features. 

In Fig. the surface density is shown along with the 
velocity field for the three CiRO-models. As a reference, 
the Roche lobe (of the relative three-body problem) is 
over-plotted. From the upper row of this figure (Mp = 1 
Mq^) we can see the patterns of the two main spirals (left 
panel). They reach the Roche lobe, but do not enter it. In 
fact, they are replaced by two streams of material which 
start from two points (located at r = 0.94, ip = 3.12 and 
r = 1.07, if = 3.07, respectively), where the flow is nearly 
at rest with respect to the planet^. Each of them enters 
the Roche lobe, encircles the planet and hits the other 
one on the opposite side. As a result of the collision, the 
material is shocked and the flow is redirected towards the 
planet. Hence, these streams assume the form of two spi¬ 
rals, winding around the planet (right panel) for 2 7 r. That 
such smaller scale spirals are detached from the main ones 
can be inferred from the direction of the flow. Along the 
main spirals the material follows the disk rotation around 
the star, moving away from the planet. Along the small 
ones the gas orbits the planet. In fact, they represent the 
outstanding features of a circumplanetary disk. A more 
detailed description of the flow regions around and inside 
the Roche lobe, concerning a Jupiter-mass planet, can be 
found in Lubow et al. ( |1999 ). 

The case Mp = 32 Mg (Fig. |^, middle row) has many 
analogies to the previous one. This planet is able to open 
a gap in the disk, as a permanent feature. However, it is 
neither so wide (the base width is 0.15 against 0.4) nor so 
deep (40% against 0.5% of the maximum surface density) 
as it is for a Jupiter-mass planet. The overall behavior of 
the matter entering the Roche lobe is very similar (left 
panel). The up-stream disk material, relative to the near¬ 
est X-point, reaches it, inverts partially the direction of 
its motion and flows into the Roche lobe. The gas stream 
penetrating from the left X-point turns about the planet, 
a,t ip < Pp, and collides with the stream incoming from 
the other X-point, generating the upper spiral arm (at 
p > Pp). However, here the locations, from which these 
gas streams depart (r = 0.97, p = 3.14 and r = 1.03, 
p = 3.135, respectively), lie closer to the LI and L2 points. 
The circumplanetary spirals are less twisted around the 


They are designated as “X-points” by Lubow et al. (1999) 



















10 


Gennaro D’Angelo et al.: Nested-grid calculations of disk-planet interaction 


1.4E-03 3.4E-01 6.9E-01 1.0E+00 9.8E-03 3.6E-01 7.IE-01 1.1E+00 


3.20 


3.15 H 


3.10- 


3.05- 



9.5E-03 


3.25 


2.4E-02 1.0E+00 2.0E + 00 3.0E+00 1.6E-01 1.1E + 00 2.1E + 00 3.0E + 00 


3.155 


3.1 60 


3.150 


3.140 


3.130 


3.120 



3.150 


3.145 


3.140 


3.135 


3.130 



0.98 0.99 1.00 1.01 1.02 


0.990 0.995 1.000 1.005 1.010 


Fig. 12. Overview of the surface density on the whole grid system for model Gino 3 (Mp = 0.5 Mq^), after 210 orbits. 
The over-plotted curve represents the Roche lobe of the restricted three-body problem. Plus signs indicate the positions 
of LI and L2 Lagrangian points, respectively on the left and the right of the planet. The color scale is linear. In order 
to enhance the contrast, in the two bottom panels S is cut at a value equal to 3 (saturated area). 
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planet than before. They wrap around it for an angle tt 
( right panel). 


For the less massive planet, Mp = 3.2 Mg, the situ¬ 
ation is somewhat different (bottom row of Fig. 13). In 
fact, within the Roche lobe, the signs of spiral fronts are 
very feeble, though some traces can still be seen. They 
assume the shape of a bar-like structure which extends, 
for 0.3 i?H, from side to side of the planet at ~ (/5p. As 
indicated by the velocity field, the circumplanetary disk 
roughly occupies the entire Roche lobe (right panel). 


Taking into account the other models as well, the fol¬ 
lowing scenario can be proposed: the lower the value of g, 
the shorter and straighter the circumplanetary spirals be¬ 
come. For example, in model Pepp 2, they track a tilde-like 
pattern, extending for a total length of 0.4 i?H- 


5.2. Density in the planet’s environment 

So far we have described qualitatively the main hydro- 
dynamic structures which are present, near the planet, 
on the length scale of a Hill radius. We shall now dis¬ 
cuss, more quantitatively, two aspects of the surface den¬ 
sity, observable on shorter length scales (< 0.5 Rh)- Both 
of them have repercussions on the torque exerted on the 
planet by the close-by matter, as explained later. First of 
all we should notice that S(r, ip) is not perfectly symmet¬ 
ric with respect to the planet, in none of our reference 
models. Such property can be checked by means of an ac¬ 
curate examination of Fig. |^. In this figure the contour 
lines of the surface density (i.e. lines on which E is con¬ 
stant) are plotted, for the finest grid levels I = ng. In 
case of CiROl (top panel. Fig. 0 , one can see that the 
right-hand arm ((p > pp) twists towards the planet more 
than the other one does. The curvature^ ratio of the more 
external spiral parts (indicated by thick lines in the fig¬ 
ure) is 0.9. Furthermore, the right-hand arm lies closer 
to the planet than the left-hand one, as can be evaluated 
from the positions of the arc centers (marked with plus 
signs). Contour line asymmetries are fainter in the case 
of Ciro2 (middle panel. Fig. 0)- However, a quantita¬ 
tive analysis shows that the external ridges of the spirals 
(thick lines) have slightly different curvatures. Moreover, 
their centers (plus signs) are not aligned with the planet 
position (represented by a big dot). In particular, if we 
consider the contours between 1.03 and 2.02 (arbitrary 
units), the right-hand ridge is a little straighter than the 
left-hand one. In fact, the ratio between their curvatures 
is 0.93. Also in this case, the arc centers indicate that ex¬ 
ternal parts of the right-hand arm are a little nearer to 
the planet than those of the left-hand arm. Indeed, sur¬ 
face density asymmetry is evident, in the case of CiROS 
(bottom panel. Fig. 0- At r ~ 1.001 and ip < ipp, the 
density is systematically lower than it is on the opposite 
side (r ~ 0.999, (p > Pp). Following the shape of density 
contours, it is possible to track what remains of the spi- 


^ The curvature of a circle of radius R is 1/R. 


ral arms. The pre-shock material conveys the more convex 
form to these lines, at r ~ 0.999 and r ~ 1.001 


5.2.1. The core 

Finally, we would like to mention what happens, in our 
models, at very shor t distances from the planet. We de¬ 
scribed, in Sect. 4.2, how material is removed from the 


neighborhood of the planet. This is usually a small frac¬ 
tion of the available matter, during each integration time 
step of the main grid. If AM®'' is the evacuated amount 
of mass, then typically 


AMe, 

M 


AE 


10 


-2 


(see the caption of Fig. |^. This choice accounts for the 
fact that a planet should not be able to accept, very 
rapidly, all the material the surrounding environment can 
offer (Wuchterl 1993 ). 

Since not all of the matter is taken away, it should pile 
up at the location of the planet, eventually forming a very 
dense core. Indeed, this is what we find, as already visible 
in Fig. 0: where density contour lines crowd around the 
planet at (cp, pp). Figure displays better how this fea¬ 
ture looks like for models CiROl and CiR02. In order to 
make a comparison of the linear extension of such cores 
with the Hill radius, we introduce the Hill coordinates. 
These are defined as a Cartesian reference frame with ori¬ 
gin on the planet and coordinates normalized to i?H (see 
the caption of Fig. for details). In these two cases, the 
core width, 2 r] (exactly defined below), can be estimated 
to be approximatively 0.1 Rh- 

One reason for the sharpness of these peaks is the very 
small length sca le w e adopt to smooth the potential of 
the planet (Sect. M)- On the finest grid level the smooth¬ 
ing length S is equal to X{ng) ~ I.4Ar(ng) (see Eq. |l^ . 
From Table one can deduce that S ~ 10“^ Rh (though 
it changes a little for the different values of q). Despite 
that, the core width always equals at least 6 Ar{ng). Other 
two hints suggest that such features are not just a nu¬ 
meric product. According to models CiR02 and WproI 
(which differ only in the number of grid levels), the core 
width is very similar. If we cut the peak at E = 2 on 
the finest level of each model, its extent, at r = rp, is 6 
and 11 grid cells, respectively. Models CiROl and Elen2 
agree in a very similar way. Eurthermore, the structure of 
these peaks does not depend on the exact placement of 
the planet within a grid cell, as models CiR03 and OiNOl 
prove. 

However, in order to study the exact properties of such 
features a more detailed physical treatment is required in 
proximity of the planet. This may be part of future com¬ 
putations which include a more detailed treatment of the 
internal constitution of the protoplanet. Anyway, proper¬ 
ties such as the local surface density profile and the veloc¬ 
ity field indicate that the core structure could be approx¬ 
imated to that of a rotating and isothermal gas, in hydro¬ 
static equilibrium. In model PeppI (1 Mj), the rotation 
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Fig. 13. Surface density and velocity field for three selected models. Top row. CiROl: grid level I = 3 (left panel) 
and I = 4 (right panel) at t = 375 orbits. Middle row. CiR02: grid level I = 4 (left panel) and I = 5 (right panel) at 
t = 325 orbits. Bottom row. CiR03: grid level I = 5 (left panel) and I = 6 (right panel) at t = 225 orbits. To avoid 
too much confusion, only 40% of the velocity field vectors are drawn. The over-plotted curve and plus signs have the 
























































Gennaro D’Angelo et al.: Nested-grid calculations of disk-planet interaction 


13 


3.160 


3.150 


3.140 


3.130 


3.120 



3.150 


3.145 

■©- 

3.140 


3.135 



0.970 0.980 0.990 1.000 1.010 1.020 


3.1 55 


3.130 


0.985 0.990 0.995 1.000 1.005 1.010 

3.148 

3.146 

3.144 

3.142 

3.140 

3.138 

3.136 


0.994 0.996 0.99B 1.000 1.002 1.004 1.006 

r 


Fig. 14. Surface density contours on the finest grid level. 
The evolutionary time is the same as in Fig. |^. Top: 
CiROl. Middle: CiR02. Bottom: CiROS. See the text for 
an explanation of the thick bow-lines and the plus signs. 




Fig. 15. Surface density plot of the highest grid level. 
Top: CiROl at t = 375 orbits. Bottom: CiR02 at t = 325 
orbits. The Hill coordinates, (xhiJ/h), are so defined: 
XB. = [j'p cos((pp) - r cos(v3)]/i?H, ?/H = [j'p sin((pp) - 
r sin(<p)]/i?H- The star is situated along the negative xr- 
axis. The azimuth (f increases with yn- 


of the core is very slow. Hence, its density profile can be 
well reproduced by that of an isothermal and hydrostatic 
gas, as demonstrated in Fig. where we compare the nu¬ 
merical results with an analytic solution for a hydrostatic 
isothermal core (solid line). 

As pointed out, the material within the core region is 
strongly coupled to the planet, due to the small distances. 
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Xh 

Fig. 16. The core width is generally much larger than that 
proper for an isothermal, hydrostatic, configuration. In 
fact the centrifugal force, due to axial rotation of the gas 
around the planet, may play an important role in sup¬ 
porting the structure. This is not the case for the core 
around an Earth-mass planet (model PeppI). Thus, the 
structure equation reads: (i<i>p = —Cg dS/E, where d’p is 
the smoothed potential of the planet. The constancy of 
the sound speed (cs = h yjG Mi^ /a) implies that the gas 
is isothermal. In the figure, the surface density of PeppI, 
at (^ = (Pp (t/h = 0), is represented by “x” signs. The 
over-plotted curve is the solution of the previous struc¬ 
ture equation. 


In some way, it may be considered as part of the pro¬ 
toplanet itself, whose structure we may not resolve well 
enough in the present paper. Whatever its nature, it is 
very likely that the angular momentum transferred, by the 
core material to the planet, may influence the planet’s spin 
angular momentum rather than its orbital one. As we are 
treating the planet as a point mass we cannot estimate its 
spin. Therefore we decided to exclude this region from the 
torque computation. To do that, we need a quantitative 
estimate of the core radius 77 , for every model in Table |^. 
We adopt the following procedure: the average density E 
is computed over the region 0.2 i?H < I’’ — J'pl < 0.5 i?H- 
Then we define the core width {2r]) to be that where 
E = 2E. In Fig. 17 (left panel) the dependence of 2 77 ver¬ 
sus q is displayed. The ratio decreases for increasing 

values of q. However, between 3.2 and 32 Mj, it seems to 
vary very little. Our measure of the core sizes is performed 
at a particular time (for CiRO-models, they are indicated 
in Fig. 1^). Anyway, such estimates are not affected much 
by this choice because the cores reach a steady state, early 
during the system evolution. As an example. Fig. [T^ (right 
panel) shows how the core mass M{r]) assumes a static 
value very soon. 


As already mentioned, the spiral features vanish for 
low values of g, and for Earth-mass planets the core be¬ 
comes the most prominent feature within the Roche lobe. 

5.3. Torque exerted on the planet 

Any protoplanet embedded in a protostellar environment 
suffers gravitational torques, exerted by the surrounding 
disk material. If the density distribution were symmetric 
with respect to the planet, or with respect to the line con¬ 
necting the star with the planet, the resultant total torque 
would be zero. However, we have just seen that this is not 
the case, and even around a planet as small as 3.2 Mj, 
the matter distribution is not axially symmetric. Thus, 
we expect a non-zero total torque from the disk, acting 
on the planet. In response to it, because of the conser¬ 
vation of the orbital angular momentum, the planet has 
to adjust its semi-major axis, which leads to a migration 
phenomenon. 

In the present computations, we evaluate the torque 
exerted on the planet in a straightforward way. First the 
gravitational force acting on the planet /g(*,i), due to 
each grid cell (i, j), is calculated. The torque, with respect 
to the star, is then 

Khj) =rpX ( 20 ) 

Since we are interested in the z-component of this vector, 
we perform the scalar product z-t{i,j) = tz{hj)j where z 
is the unit vector of the vertical direction. Finally, by sum¬ 
ming over all i and j, we obtain the total disk torque 7 d, 
whereas the simple contraction over the azimuthal index 
j gives the radial distribution of the torque. 

The quantity tz{i,j) is computed on each grid level. 
Where the computational domain is covered by more than 
one level, the torque tz{i,j) on. the finest grid is considered 
for the evaluation of 7 d . 

In order to avoid the region dominated by the core, 
we exclude a certain area from the computation of Tu- 
Because of the way we defined the core radius 77 , some 
core material may still lie outside |r — rp| = 77 . Therefore, 
for safety reasons, we choose not to take into account the 
planet neighborhood defined by |r — rp| < /3 = 2 77 . The 
only level, upon which this operation is relevant, is the 
highest since it provides the gravitational torque form the 
regions closest to the planet. On coarser levels such op¬ 
eration is meaningless since inner torques are taken from 
elsewhere. But it can be useful to confer a more regular 
look to the radial torque profile. We generally adopt the 
value /? = 5 6 {l) for I < ng. 

However, as this might be somewhat arbitrary, we dis¬ 
cuss in a separate section how the choice of /3, on the finest 
level, affects the value of the total torque Td- 

Because of the global angular momentum transfer, the 
disk material (orbiting the star), at r > rp, exerts a neg¬ 
ative torque on the planet whereas the inside-orbit gas 
tends to increase its angular momentum. This tendency 
changes, as material closer to the planet is accounted for. 
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Fig. 17. Left panel. Core width, 2r/, in units of the Hill radius as a function of the mass ratio q. CiRO and Pepp- 
models are considered along with GiN02 and GiNOS. The measure errors are assumed to be equal to ^ (see Table p. 
Righ panel. Mass within a distance rj from the planet versus time. The core mass M(rj) is normalized to the mass of 
the planet. 


and it may reverse eventually, once in the circumplanetary 
disk (for Mp = 1 such behavior was also found by 

Lubow et al. 1999 ). The radial torque distribution, from 
GiRO-models, is shown in Fig. H- The two sets of pro¬ 
files belong to the grid levels I = ng — 1 (top row) and 
I = ng (middle row). In the case of GiROl (top row, left 
panel), the sign reversal of the torque is not completed 
yet. However some negative torques are exerted from re¬ 
gions inside the planet orbit and some positive torques 
are exerted from the opposite side. On the domain cov¬ 
ered by this grid (fourth level), the torque contribution 
coming from cch < 0 is positive while the one coming from 
a;H > 0 is negative. The resulting net torque is positive, 
as the magnitude of the latter contribution is 2.4 times 
smaller than that of the former. The torque behavior gets 
more complex if we restrict to a region closer to the planet 
(middle row, left panel). Though not evident at a first 
glance, the signs are reversed if compared to the preced¬ 
ing grid level. The torque exerted by the region xh > 0 
is definitely positive and, in magnitude, almost 30 times 
as large as that arising from the region xh < 0. Thus, 
this region exerts a strong, positive, net torque. Indeed, 
the phenomenon of the torque sign reversal is clear in the 
case of Giro2 and CiROS (Fig. [1, center and right pan¬ 
els). For both grid levels, inside-orbit material lowers the 
angular momentum of the planet while outside-orbit ma¬ 
terial acts in the opposite direction. On the finest level 
of Giro2, the ratio of the negative to the positive torque 
contribution is just 0.96 (in absolute value), whereas it is 
0.3 for Giro3. 



i 


Fig. 19. Partial torque T(/), normalized to |7 d|, plotted 
against the grid level I for: GiROl, GiR02 and GiR03. 
T{ng)/\TY^\ = —1 because, by definition, T{ng) = 7d and 
7 d is negative for these models. 


In order to check which is the overwhelming contribu¬ 
tion, between positive and negative torques, on the vari¬ 
ous grid levels, we can define the partial torque Til). This 
represents the total torque computed over the entire do¬ 
main but the part covered by the grid level I + 1 and such 
that ^{ng) = 7d. Fig. illustrates the sign of the par- 
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Fig. 18. Radial distribution and two-dimensional contour map of the gravitational torque exerted by the disk on the 
planet. Left panels. CiROl: levels I = 4 and 1 = 5. Center panels. CiR02: levels I = 5 and 1 = 6. Right panels. 
Ciro3: levels I = 6 and I = 7. The two-dimensional torque distribution (bottom panels) is reported for the finest grid 
level and is normalized to its maximum, absolute, value. 


tial torque and its relative strength, on each grid level, for 
CiRO-models. We can see that the total torque is negative 
in all of the three models. In CiROl and CiR02 all levels, 
but the highest, contribute to lower the planet angular 
momentum. On the contrary, the matter inside the finest 
level raises the overall torque (of a fair amount in case 
of CiROl). In Ciro 3, a positive torque is exerted by the 
material outside a region, around the planet, with a linear 
extension ~ 2hrp. Instead, levels 3, 4, 5 and 6 provide 
negative torques, which are then weakened considerably 
by the positive torque coming from the finest level. 

An overview of the torques exerted by different por¬ 
tions of the disk, for all the relevant models, is given in 
Table |[ First of all we notice that, within 0.5 i?H, the 


phenomenon of the sign reversal is observed in all of the 
models, but PeppI and GiN02. Anyway, in the latter case, 
the torque inside this region is positive. 


In CiRO-models and in model GiN03, the reversal of 
the sign is such to produce a positive torque over the do¬ 
main |r —Tpl < 0.5 i?H- This is because the surface density 
of the leading material is slightly higher than that of the 
trailing matter (Sect. 5.1 and Fig. [T^ ). In contrast, the 
torque is negative in case of Pepp2, Pepp 3 and Pepp 4. 
We note, however, that the positive torques exerted by the 
outside-orbit material, within this region, strongly atten¬ 
uate the magnitude of the negative net torque exerted by 
the rest of the disk. 
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Table 2. Gravitational torque exerted on the planet aris¬ 
ing from different disk regions, for the relevant models 
of Table ^ The entire domain is divided into three re¬ 
gion: outside the Hill circle; inside the circle of radius 
i?H/2 and the zone in between. Then they are divided 
further in order to distinguish between inside-orbit {in) 
and outside-orbit {out) contribution. The torque is com¬ 
puted over each of these regions, by employing (3 = 2 rj and 
according to the value of the normalized distance from the 
planet d =\r — rp|/i?H. Torques are normalized to 7 d. 


Model . 

d < 0.5 

0.5 < d < 1 

d > 1 

in 

out 

in 

out 

in 

out 

CiROl 

-0.03 

2.08 

-1.18 

-0.43 

0.53 

-1.97 

Ciro2 

-3.07 

3.23 

-1.06 

0.94 

2.04 

-3.08 

Ciro3 

-1.93 

3.13 

-0.49 

-0.97 

7.26 

-8.00 

PeppI 

0.33 

-0.40 

0.33 

-0.36 

4.22 

-5.12 

Pepp2 

-2.32 

2.05 

-0.18 

-0.18 

1.91 

-2.28 

Pepp3 

-2.61 

2.07 

-0.16 

-0.01 

0.85 

-1.14 

Pepp4 

-3.01 

2.41 

-0.29 

0.25 

0.84 

-1.20 

Gino2 

2.47 

-1.65 

-1.44 

1.02 

2.56 

-3.96 

Gino3 

-1.54 

5.12 

-2.06 

-0.34 

1.90 

-4.08 


Clearly, since neighboring material may tend to reduce 
the magnitude of negative torques acting on the planet, 
it could either slow down its inward migration or reverse 
the direction of its motion. 

As anticipated, the sign reversal of the radial torque 
distribution is not observed in model PeppI (Mp = 1 
Mg). Due to the very low mass of the perturber, the only 
structure present inside the Roche lobe is the density core. 
Whatever level is considered, the inside-orbit gas always 
provides positive torques while negative ones come from 
the outside (see also Table Negative torques are some¬ 
what stronger, on any grid level. Almost the 50% of the 
total torque is generated between ~ O.Sfirp and ~ 
at the starting positions the disk spirals. 


5.3.1. 2D-torque distribution 

The detailed balance of the torques, arising from very close 
material, depends on the medium and small-scale density 
structures around the planet, such as the shape of spirals. 
Therefore, referring to what was stated in the previous 
section, we can deduce that the radial torque asymmetry 
is a direct consequence of the asymmetric distribution of 
the gas with respect to the planet. 

Thus, a more comprehensive description of the torque 
behavior requires its full two-dimensional distribution. In 
the bottom panels of Fig. |^, the contour lines of the 
two-dimensional torque, tz, are shown for each reference 
model. The interesting point to notice here is that the 
largest magnitude torques arise from the corotation loca¬ 
tions, i.e. where r = r-p (a;H — 0). Here in fact, fg{i, j) 
is perpendicular to Vp and the cross-product in Eq. ( P0[) 
achieves its maximum (minimum). The material leading 
the planet (at (p > (pp or j/h > 0 ), pulls it ahead and 


makes it gain angular momentum. The trailing material 
brakes the planet making it lose angular momentum. 

Let’s consider two fluid elements at (0, ± yn) and write 
their mass density as Then we can write |tj| oc 
TpH^/y'^, which yields: 


, , AE 

IG I rp —, 

Vh 


( 21 ) 


where AE = E+ — E“. Any mismatch of the surface 
density, AE, causes a torque mismatch amplified by an 
amount equal to It’s worth noticing that, on larger 


distances \r — Vp 
less, in fact: 


|t+| - IG I ot 


AE 


Tp, the torque mismatch is amplified 


( 22 ) 


This is the reason why surface density asymmetries near 
the planet have a very strong impact on Tjj, and they can 
easily prevail against the more distant ones. 

The region responsible for the maxima and minima of 
the radial torque distributions in the middle row of Fig, [l^ 
can be identified by means of the 2D-torque maps. For ex¬ 
ample, in the case of CiROl, we see that the maximum at 
XH = —0.15 and the minimum at xh — 0.15 are produced 
at yn — 0.1 and j/h — —0 .1, respectively. In the other two 
cases, radial distribution extremes rise from regions where 
the torque function tz is steeper than it is on the opposite 
side of the planet. 


5.4. Planet migration 

If we consider a planet moving on a circular orbit, we find 
that the rate of change of its semi-major axis a, caused by 
an external torque Tq, is 


da 

dt 


2Tr, 


Mp a Op 


(23) 


Analytical estimates of Tq show that two limiting cases 
exist, depending on whether the planet is massive enough 
to generate a gap or not. In the first case of small plane¬ 
tary masses, we have the so-called type I migration. Ward 


(1997) derives the following expression: 


da\ 1 , ■? 

— --g h~ ^ 

dt Jp 2 


O,, 


M4 


(24) 


The direction of the migration is inwards because of the 
dominating role of the outer Lindblad resonancesn (Ward 


1986, 1997). In the second case the planet is more massive, 
opens up a gap, and the evolution is locked to that of the 
disk. As a general trend, the disk material drifts inwards 
on the viscous time scale, and so does the planet. It follows 
that 


da \ 3 3,0^ 

I =-= — ah a PL- 


dt /j 


2 a 


P’ 


(25) 


This is true as long as th e temperature gradient within the 


disk is negative (Ward 1986). 
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Fig. 20. Migration time scale tm versus the mass ratio q. 
Open triangles represent the results from models Giro, 
Pepp, Gino2 and Gino3. The total disk torque 7d is 
computed assuming /3 = 2 77 , for each model. This means 
that the region lying inside |r — rp| = 2 rj is not taken 
into account. In order to express tm into physical units, 
we suppose that the planet orbits at a = 5.12 AU in a 
disk with mass = 3.5 x 10“^ The curve over¬ 
plotted represents predictions of the analytical theory as 
formulated by Ward (1997), for the case of “strong” vis¬ 
cosity (a = 4 X 10“^ in our case) and accounting only for 
Lindblad torques. It is derived assuming an unperturbed 
constant surface density and a disk temperature dropping 
as 1/r. The behavior of this curve reduces to Eq. ( p^ 
and Eq. (p5|) letting g —> 0 and g —> 00 , respectively. The 
vertical line marks the value given by Eq. (pq). 


In Fig. p^the migration time scale tm is shown for the 
main models listed in Table The drifting motion, is di¬ 
rected inwards, in all cases. The lowest migration velocity 
belongs to the Earth-mass planet (PeppI). The second 
lowest drifting velocity is that of the Jupiter-mass planet 
(GiroI). The most rapidly migrating planet is the one 
having Mp = 20 Mj (Pepp4). 

In agreement with predictions of analytical linear the¬ 
ories, the drift velocity |a| increases for increasing planet 
mass, just as prescribed by type I migration (Eq. The 
fast speed branch has a turning point around g ~ 6 x 10“^, 
after which migration slows down considerably. Past this 
point, |a| drops as the planet mass increases. According 
to the linear theory, this property announces the transi¬ 
tion to type II migration (Eq. p^. As a comparison, the 
complete theoretical behavior of tm is also reported in 
Fig. p0|. It was derived by War d (|I997 ) for viscous disks 
with a > lO”"^. Eqs. (p^ and (|25|) represent the asymp¬ 
totic branches of this curve for very light and very heavy 
planets, respectively. The vertical line, in the figure, indi¬ 
cates the mass value above which migration is faster than 
inward viscous diffusion. 


In all the models under examination, numerical simu¬ 
lations predict a slower drift than analytic linear theory 
does. Roughly, migration is two times as slow for models 
Gino2, Gino3, Pepp3 and Pepp4; three times as slow 
for GiroI, Giro2, Giro3 and Pepp2. Planets as massive 
as 3.2 and 64 Earth-masses migrate on the viscous time 
scale of the disk (Eq. p^. tm is six times as large in case 
of the Earth-mass planet, when compared to analytic pre¬ 
dictions. Such discrepancies are likely to arise also because 
the theoretical curve in Fig. does not include corota¬ 
tion torques. Yet, we have seen that just these torques 
slow down inward migration, preventing the total torque 
7 d from being very negative. 


It’s worthy to note here that the relevance of coorbital 
torques, to the orbital evolution of a protoplanet, was al¬ 
ready pointed out by Ward (1993). 


which is known as type II migration. Gomparing Eq. (p^ 
with Eq. (p5|), it turns out that type I drift is faster than 
type II (i.e. faster than viscous diffusivity) whenever 

The parameter values employed herep] yield a right-hand 
side equal to 2.8 x 10“®, which is just a bit smaller than 
an Earth mass (0.93 M^). Fast type I migration should 
continue till the planet grows enough to impose a gap on 
the disk. By that time, however, it could have already 
impacted the parent star. Once entered this fast drifting 
regime, it seems that the planet may survive only if the 
growth time scale tq = Mp/Mp is much smaller than the 
migration time scale tm = a/\a\. 

® If we take into account the dependence of the unperturbed 
surface density upon the radial distance r, Trad'S — 0.38Md. 


5.4.1. /3-dependence 


In the above discussion, as explained in Sect. 5.3, we do 
not account for the torque exerted by matter closer than 
/3 = 277 to the planet. Now we would like to relax this 
assumption and consider how material, lying even closer, 
affects the overall torque 7h- We mentioned already that 
the tendency of the nearby gas is to increase the angular 
momentum of the planet. As we enter the core dominated 
zone, such tendency may grow stronger and stronger, re¬ 
ducing more and more the magnitude of the negative 
torque exerted by the rest of the disk. 

An overview of the effects, due to nearby matter, on 
the migration time scale tm, is given in Fig. Actually, 
these plots show the dependence of ti/tm oc Tp upon /3, 
where v = d/|d| indicates the direction of the planet’s 
migration. The distance where /3 = 2 77 is marked with 
a vertical, solid line. 7 d is directly proportional to the 
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Fig. 21. Migration time tm as a function of /3, the radius 
of the region excluded from the computation of the torque. 
The ratio v = a/\a\ indicates the direction of the planet 
drift: u < 0 for inward migration. The semi-major axis a 
has the same value as in Fig. ^ The disk mass is cast into 
the form Md = k x 10“^ Tf*. Since in this paper we are 
using the value k = 3.5, the factors in round brackets be¬ 
come 2.9 X 10^ and 2.9 x 10^ years, respectively. The solid, 
vertical, line indicates the length /3 = 2 77 (see Fig. 0 )- 


mass of the disk Md. Therefore, to remove this potential 
restriction, in these plots we let it as a free parameter and 
write Md as k x 10“^ M*, although we use k = 3.5 for 
our estimates. In Fig. ^ 7 d = ThiP) is presented for all 
of the three reference models. 

In CiROl (top panel), 7 d becomes larger as (3 gets 
smaller. The sign of the total torque changes around 
(3 = T], the threshold of the density core. As a compar¬ 
ison, the behavior of 7 d, versus (3, is reported also for the 
model ElenI. In this case, because of the closed inner ra¬ 
dial border, the amount of matter inside the orbit of the 
planet is five times as large as that of CiROl. Outside the 
Hill circle, the torque exerted by the inner-disk, in case 
of ElenI, is twice as large as that measured in CiROl. 
Instead, torques arising from the outer-disk nearly coin¬ 
cide. Inside the Hill circle, in ElenI, the contribution to 
7 d is relatively small down to ~ 0.1 whereas, in CiROl, 
it never appears to be negligible. 

Ciro2 (middle panel) behaves somewhat differently 
from CiROl. The total torque attains a minimum around 
(3 = 0.15 i?H, where tm — 3 x 10"^ years. Then the pos¬ 
itive torques, exerted by close matter, increase the total 
torque, though it remains negative all the way down to 
(3 = 0.03 i?H. Below such value, 7 d diverges positively. 
Results from the higher resolution model, WproI, do not 
differ significantly (dash-line in Fig. ^ ). 

7d varies smoothly, as a function of /3, in case of 
Ciro3. The torques arising from the region enclosed be¬ 
tween f3 ~ 0.2 i?H and (3 « 0.4 i?H almost cancel out, so 
that the total torque appears nearly constant (tm « 10 ® 
years). At shorter distances, positive torques prevail over 
the negative ones and 7 d starts to increase. The sign of 
the total torque reverses at /3 ~ 0.12 Ru. For example, 
at P = 7], its value is quite positive, imposing an outward 
migration rate tm — 3 x 10^ years. 

Some comments should be devoted to how the smooth¬ 
ing length affects the total torque. We did not try to re¬ 
duce further its value, however we ran some models, iden¬ 
tical to CiROl, but multiplying A (see Eq. ^ by some 
integer number. A larger smoothing length tends to smear 
out more the surface density nearby the planet. Besides, it 
also causes the material to be distributed more symmetri¬ 
cally around it. Both tendencies contribute to reduce the 
magnitude of the net torque arising from a region with 
radius w A. 

5.5. Circumplanetary disk: gas flow 

In Sect. H we have described in detail how gas flows into 
the Roche lobe of a planet, for three particular values of 
its mass. 

If a planet is massive enough, say Mp > 10 Mg, 
the streams of matter, entering the Roche lobe, produce 
strong shock waves which then rule the gas flow inside 
this region. Material passing through the shock fronts is 
deflected towards the planet, tightening its orbit on it. 
Less massive planets are not able to cause strong pertur- 
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Fig. 22. Contour lines of velocity ratios rCm/wS (top panels) and Wrot/Wrot (bottom panels). From left to right, we 
report the above quantities, as computed on the highest level, of models: CiROl, CiR02 and CiROS. Win > 0 contours 
indicate places where the flow approaches the planet. At locations where the in-fall velocity is zero the flow performs 
a pure rotation around the planet. 


bations inside the Roche lobe. As a consequence, the flow 
pattern appears more uniform around the planet. 

Now we would like to investigate quantitatively the 
rotational regime of the gas inside the circumplanetary 
disks. In particular, we would like to estimate how much 
it resembles a Keplerian one. In order to address this issue 
we decompose the local velocity field u, in two compo¬ 
nents, representing the in-fall velocity Win and the rota¬ 
tional velocity iCrot of the fluid relative to the planet. The 
first component is defined as: 


Win = -U ■ 


which is positive if a fluid element moves towards the 
planet. The quantity Wrot is the projection of u along the 
direction orthogonal to r — r p and such that (lUrot x Win) • ^ 
is positive. With this choice Wrot is positive for a counter¬ 
clockwise rotation. 


If circumplanetary disks were regular accretion disks, 
we should expect them to be in a “Keplerian” regime. This 
is characterized by the rotational velocity 





(27) 


and the inward viscous diffusion 

D ^ ^ 

w- = 


2 r - r. 


(28) 


Equations (|^) and ( |28| ) don’t include the smoothing 
length, (5, because its effects were checked to be unimpor¬ 
tant. Comparing Wmt and wu 


with Eq. (27) and Eq. 
respectively, we can estimate how much the circumplan¬ 
etary flow is close to be Keplerian, i.e. close that of an 
unperturbed viscous disk. 

Figure ^ shows, for CiRO-models, the contour lines of 
Win normalized to wj^ (top panels) and w^t normalized to 
(bottom panels). 

As first remark we note that, if we compare Fig. to 
Fig. 1^ lines of equal surface density perturbation are also 
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Orbits 


Fig. 23. The plot shows the ratio Wrot/w^^ at (p = (pp 
{uh = 0 ), as computed on the grid level ^ = 6 of the model 
Elen2 (Mp = 1 Mqp). The core area, in the Jupiter-mass 
case, extends for [0.08 x 0.08] Rb^ (see Fig. left panel). 
The resolution, in Elen2, is such that this is covered by 
~ 12 X 12 grid cells. 


Fig. 24. Mass accretion rate onto the planet versus time, 
in units of 10 “'^ disk masses per orbital period of the 
planet. For Md = 3.5 x 10“^ M*. and a = 5.2 AU, one 
dimensionless unit corresponds to 2.95 x 10“^ 

The initial accretion rate is very small for model CiROl 
because of the imposed initial gap. 


lines of equal velocity perturbation, as spiral wave theory 
predicts. 

From the top panels of Fig. ^ we can see that ma¬ 
terial approaches the planet along well defined patterns. 
Contours u>in = 0 mark locations where the flow rotates 
around the planet without altering its distance from it. 
They also separate regions in which material proceeds to¬ 
wards the planet from those where it moves away. One of 
these contours runs along the spiral ridges. Across it, the 
in-fall velocity changes abruptly its sign. 

The ratio ?Cin/w£ becomes smaller as the gas comes 
closer to the planet. Since the viscous diffusion wP is not 
related to Mp, it’s possible to compare the magnitude of 
Win for the different cases. Contour level values indicate 
that it gradually reduces as Mp gets smaller. 

As regards the rotational component of the veloc¬ 
ity field Wrot (Fig- bottom panels), we can see that 
it is generally slightly below For the Jupiter-mass 
case, this can be seen also in Fig. Centrifugal over¬ 
balance regions are not present around the smallest planet 
(Mp = 3.2 Mg). Instead, they are observed in CiROl, at 
(xhjI/h) « (—0.1,0.3) and (0.15, 0.1), and in CiR02, where 
they are labeled. In both cases, anyway, the ratio Wrot/'R'rot 
is very close to unity. Centrifugal under-balance is mainly 
established along spiral ridges. This is in agreement with 
the idea that spiral arms are zones of compression, hence 
pressure plays a more active role in supporting the gas. 


In all of the three cases shown in Fig. ^ w^ot reveals 
somewhat a circular symmetry only within a distance ~ 77 
from the planet. Yet, we have to notice that this coincides 
nearly with the region from which matter is removed to 
simulate the gas accreti on. F igure also indicates that 
the core material (Sect. ^.2.1 ) rotates slower, when it ap¬ 
proaches its center. 


5.6. Accretion onto the planet 

Gas matter closely orbiting the planet is eligible to be 
accreted once its distance, |r — rpj, is less than Kac « 
9 X 10“^ i?H. The details of this process are described in 


Sect. O (see Fig. In general, the mass accretion rate 
of a planet, Mp, becomes relatively constant after the gap 
(if any) has evolved to a quasi-stationary state. 

For a Jupiter-mass planet (CiROl) this happens 
around 100 orbital periods, as indicated by the solid line 
in Fig. The theoretical gap imposed at the beginning of 
the evolution (see Fig. 0 ) is deeper and wider, at p = pp, 
than it is later on. For this reason Mp is negligibly small 
at early evolutionary times. The partial replenishment is 
related to the formation of the circumplanetary disk which 
supplies matter for the accretion process. 

Smaller planets dig narrower and shallower gaps so 
this quasi-steady regime is reached even earlier. For both 
Ciro2 (Fig. 1^, short-dash line) and CiR03 (long-dash 
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line), Mp reduces a little as the evolution proceeds. This 
is likely due to the depletion of the inside-orbit disk. In 
case of CiROl, most of the inner-orbit material is cleared 
out during the transitional phase (75% after 100 orbits). 
Therefore, it does not contribute much during the quasi¬ 
steady phase. 

In Fig. the accretion rate is plotted against the 
mass ratio (left panel), in dimensionless units as in Fig. 

Mp increases as the planet mass increases and reaches a 
maximum around 0.5 M-j^ (model Gino3). Beyond this 
point the curve starts to drop. In fact, the accretion rate 
of CiROl lies between those of Ciro2 (Mp = 32 Mj) 
and Gino2 (Mp = 64 Mg). The accretion rates for plan¬ 
ets with masses above about 10“'^ M* (~ 30 Mg) are 
consistent with those obtained by recent models studying 
the evolution of giant proto-planets (Tajima & Nakagawa 


1997). For smaller masses our accretion rates are substan¬ 


tially higher. However, all the more detailed studies of 
protoplanetary evolution are spherically symmetric (see 
the review by Wuchterl et al. 200C), and accretion via a 
planetary accretion disk may allow for higher accretion 
rates. 

We did not perform computations involving planets 
heavier than one Jupiter-mass so we cannot follow the 
trend of the curve for larger values of q. However, Lubow 
et al. (1999) found that Mp decreases in the mass range 
from 1 to 6 Mq^. 

In the right panel of Fig. the growth time scale 
TG = Mp/Mp is plotted versus q. Logarithmic scaling of 
the axes shows that the curve decreases almost linearly 
with respect to the mass of the planet. If we perform a 
linear least-square fit of the values for which g > 2 x 10“^, 
corresponding to Mp = 64 Mg, we get: 


= - 0 . 66 . 


dl0g(l/TG) 
dlog(g) 

This equation yields: 


0.66 2/3 

TG oc g ~ g ^ . 


(29) 


At higher values of g, the curve steepens, decreasing more 
rapidly. Taking into account the growth time scales rela¬ 
tive to the two most massive planets, one finds 


TG oc gl-34 g4/3_ 


(30) 


Therefore, as a planet becomes more massive, the growth 
time scale grows with an increasing power of its mass. 
Thus, we can argue that very high mass planets should be 
extremely rare. 


5.6.1. Influence of numerics 

Finally, we comment on the influence of some numerical 
parameters upon the mass accretion rate. Table lists 
some results, after the same number of orbits, for various 
models which can be useful to the goal. The value of Mp 
may depend on either the radius Kac of the accretion region 
and the mass evacuation rate Kev By comparing models 


Table 3. Mass accretion rate onto the planet for different 
numerical parameters. This quantity is reported at the 
same evolutionary time for similar models. Mp is given in 
units of disk masses per orbit. 


Model 

g 

Mp 

Orbits 

CiROl 

1.0 X 10’^ 

8.7 X 10”® 

200 

Elen2 

1.0 X 10"® 

8.9 X 10“® 

200 

Ciro2 

1.0 X lO""* 

9.0 X 10“® 

280 

WproI 

1.0 X lO""* 

8.6 X 10“® 

280 

Wpro2 

1.0 X lO""* 

1.0 X lO-'^ 

280 

Ciro3 

1.0 X 10"® 

4.5 X 10“® 

120 

GinoI 

1.0 X 10"® 

4.5 X 10“® 

120 


Ciro2 and Wpro2, it turns out that doubling Kgv, the 
accretion rate is only 11% higher. In order to estimate 
how relevant the extension of the accretion region may 
be, we ran a model having an accretion length scale (Kac) 
0.6 times as small as that of GiN03 (the rest of the model 
set-up being identical). We obtain an accretion rate 7% 
smaller. 

Further, Mp could depend on the numerical resolu¬ 
tion as well. Indeed, this dependency turns out to be 
very weak, as indicated by two sets of models. Elen2 
and Wpro2 have, in the accretion region, a resolution 
two times as high as that of GiROl and CiR02, respec¬ 
tively (see Table 1^). Despite this fact, in the first case, the 
accretion rates differ by just 2%, whereas the difference 
amounts to 4% in the second. 

As mentioned in Sect. 4.4, the planet in not symmet¬ 
rically centered within a grid cell. Model GiNOl was de¬ 
signed to accomplish a fully symmetric configuration (as 
already explained in Sect. ^^). However, the planet ac¬ 
cretion rate is not affected by this position shift, as can 
be seen by comparing the values in Table |^, belonging to 
GinoI and CiR03. 


6. Conclusions 

A number of numerical simulations concerning disk-planet 
interactions have been performed to get new insights into 
the scenario of the joint evolution of protoplanets and their 
environment. They have confirmed analytical theories for 
gap formation and planet migration. 

However, many open questions still remain. The most 
important unsolved issue is the influence of the ambient 
gas on the dynamical evolution of a planet. Another one is 
the way disk-planet interaction changes when small plan¬ 
ets, in the mass range of Neptune and Earth, are consid¬ 
ered. 

We began to investigate in both directions by means of 
a nested-grid technique, which is particularly suitable for 
treating these problems. The main asset of this numerical 
scheme is the possibility of achieving, locally, a very high 
spatial and temporal resolution. With such a method we 
are able to resolve very accurately both, the inner parts 
of the Hill sphere of the planet and the global structure 
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Fig. 25. Left panel. Mass accretion rate Mp as function of the mass ratio q. Quantities are normalized to the disk 
mass Md. Right panel. Dependence of the growth time scale tq = Mp/Mp on the planet mass (filled triangles). 
Md = 3.5 X 10“^ M* is employed to express tq into years. The two-branch solid line represents the first-order 
polynomial fitting the data (see Eqs. and ^). 


of the disk, because we treat the whole azimuthal extent 
of the disk. 

Concerning the issues mentioned above, with the 
present paper we tackled some outstanding problems con¬ 
cerning the growth and migration of protoplanets, cover¬ 
ing a range from one Jupiter-mass down to one Earth- 
mass. Thus, even though we do not include the detailed 
energetic balance of the planetary structure, which is still 
beyond present day computer facilities, this study repre¬ 
sents a definite improvement in the determination of the 
torque balance on protoplanets. 

Our main achievements can be summarized as follows. 

1. Inside the Hill sphere of the planet a circumplanetary 
disk forms. Within it strong spiral density waves de¬ 
velop, if the planet mass is larger than ~ 5 Mg. These 
waves assume the shape of a two-arm pattern. The two 
spiral arms are slightly asymmetric with respect to the 
planet. For decreasing planet masses, they stretch and 
shorten. Matter is observed to pile up at the location 
of the planet, generating a very high density zone (we 
named as density “core”), which might represent its 
primordial gaseous envelope. 

2. Nearby material exerts positive torques on the planet, 
slowing down, considerably in some cases, its inward 
migration. Most of these torques arise from corota¬ 
tion regions, i.e. from gas lying on the planet’s orbit. 
Analytical models about migration do not account for 
them. This can be one of the reasons why our estimates 
of the migration time scales give somewhat higher val¬ 
ues than those predicted by such theories. 


3. Within a distance of ~ 0.1 from the planet, the 
point-mass approximation becomes too restrictive and 
maybe not appropriate. Therefore, the structure of 
the planet should be also taken into account over 
such a length scale. This is absolutely necessary if one 
wants to evaluate how much of the angular momentum, 
transferred by closely orbiting matter, is conveyed to 
the spin of the planet rather than to its orbital angular 
momentum. 

4. The Keplerian rotational regime of circumplanetary 
disks is affected by spiral perturbations. Just as for 
the mass density, the more massive the planet is, the 
stronger such perturbations are. Gas material, pass¬ 
ing through the spiral fronts, is deflected towards the 
planet. Instead, in the inter-arm regions it moves away 
from it. This is in analogy with the spiral wave theory 
in galaxies. Around Earth-mass planets, the rotation 
of the gas is very slow if compared to the Keplerian 
rotation. In fact, in this particular case, the density 
core has nearly a hydrostatic structure. 

5. The mass accretion rate, as a function of the mass of 
the planet, has a maximum around Mp = 0.5 Mtj,. 
As long as Mp < 0.2 Mq^, the growth time scale of a 
planet increases, approximatively, as Mp^'^^. For more 
massive planets, it increases roughly as Mp'^^^. Such 
a dependence may contribute to limit the size of a 
massive planet. 
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Since we just started to explore these new grounds, each of 
the items above may deserve a more specific and dedicated 
study. 

Next efforts should be devoted to refine the physical 
model, especially in the vicinity of the planet. Henceforth, 
some of the future developments could be: 

— including an energy equation, by implementing an ap¬ 
proximate treatment of radiative transfer and viscous 
dissipation; 

— improving the equation of state by using an alternative 
form which accounts for the planet’s structure; 

— evaluating possible effects due to the two-dimensional 
approximation of the disk, via three-dimensional sim¬ 
ulations. 
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